% This Runs the ABM multiple times, randomly adjusting model parameters,
% choosing the best model, and repeating a desired number of times or until
% the fit to historical data reaches a desired threshold

% Attempting to get best case, especially focusing population size

%% Set Up
% Settings
OnlyGoodData = 1; % Set to 1 if you want cases filter to those deemed "Good"
PlotRaw = 0; % Set to 1 to plot raw data for each Good case
LogisticOrNot = 0; % Set to 1 for Logistic Regression, 0 for Rolling Average
LogOrNot = 1; % Set to 1 for logged peak participation, 0 for not
MatchUpTo = .5; % When calculating error, only calculate up to this value for X
ncases = 25; % Number of cases to generate each loop
MaxIter = 4; % Maximum number of loops
R2PopGoal = .2; % Good Enough Population R2
R2Goal = .95; % R Squared Goal to stop looping 

%Get NAVCO 1.2 Data
HistoricalData = GetNAVCOData(4,LogisticOrNot,LogOrNot,0,[],0,0); %Nonviolent Data, Logistic, Log, No plotting
PopulationBins = [-1:.5:1.5];
temp = GetRollingAverage(HistoricalData.XRaw,HistoricalData.YRaw,PopulationBins,1,.25,1);
HistoricalData.BinNumber = temp.nsize;clear temp;

%% DeclareSimulation Parameters
% Simulation Variables
sim.PauseTime = 0; % Seconds to pause each time step
sim.PlotOn = 0; % 1 to make plots, 0 for no plots
% Parameters
clear p
p.MaxSteps = 200; % Maximum number of time steps
p.LatticeX = 40; % Number of spaces in graph in X direction
p.LatticeY = 40; % Number of spaces in graph in Y direction
p.torus = 1; % Should lattice connect at sides? 1 Yes, 0 No
p.PercentFillCitizens = 70; % Percent Lattice is filled with agent type
p.PercentFillActivists = nan; % Percent Lattice is filled with agent type
p.PercentFillPolice = 4; % Percent Lattice is filled with agent type
p.vision = 4; % Number of spaces an agent can see in each direction
p.MaxJailTerm = 10; % Max Jail Term in time steps
p.StartingGovernmentLegitimacy = [.6 .05 .1]; % Starting value (between 0 and 1) for government legitimacy
p.ChanceFindNVResistor = 40;
p.ChanceTargetNonviolent = 25; % Percent chance police will target a nonviolent agent
p.ChanceKillNonviolent = 10; % Percent chance police will attempt to kill nonviolent agent
p.BackfireCoefficient = .99; % Factor by which government legitimacy decreases
p.f = [.07 .002 .05]; % Threshold used for determinging behavior
p.ProtestCycle = 5;%[10.5 3 3];
p.ProtestDuration = 1;%[1 1 1];
p.nNV = [1 .5 1];[.02];
p.DelayStartMax = [5];
p.PeerPressureNumber = [3.46 1 1];% # of Protestors for PeerPressure to equal .5 (and vision = 7)
p.PillarProxStrategy = 0;
p.ActivistSearchVision = 10;
p.PercentFillPillars = [1.0 .3 .5]; % Percent Lattice is filled with agent type
p.DefectThresholdStdDev =0; % Defect Threshold does not vary pillar to pillar
p.DefectThreshold = [.05 .03 .05];
p.DefectThresholdMean = nan;%.05;%[.05 .002 .05]; % Threshold for pillars/police to defect
p.DefectThresholdSTD = nan;%[.02 .02 .005]; % Threshold for pillars/police to defect
p.NonviolentSuccessPercent = [1 10 1];
p.NVVictoryPercentMean = nan;1;%[1 .5 1];
p.NVVictoryPercentSTD = nan;[10 3 1];
p.PillarProxStrategyMean = nan;[.5 .05 .3];
p.ReorderAgentsParam = 1;
p.PercentFickle = [95 5 70]; % Percent of Citizens who join and rejoin resistance based on seeing protest. The nonfickle stick remain regardless;
p.PercentImmediateProtest = [20 10 0]; % %Percent of Citizens who will immediately protest after joining


% Create Matrix of All Possible Variables
[ParamCases] = param4evolution(p,ncases);

%% Start Loop
LoopCondition = 1;
LoopCounter = 1;
R2History = [];
NewParamCases = [];
while LoopCondition
    % Update ParamCases if not first iteration
    if LoopCounter > 1
        ParamCases = NewParamCases;
    end

    %% Run Cases
    %Sets Values to Be Varied each run of a given case
    PopulationMean = .8;
    PopulationStdDev = .3;
    PopulationCases = 200;
    vary = randn(PopulationCases,1)*PopulationStdDev+PopulationMean; %Varies initial activist density
    vary(vary<=.1) = .1;
    vary2 = (randn(PopulationCases,1)).^2; % Needed For Chi Squared Distribution
    vary3 = (randn(PopulationCases,1)).^2; % Needed For Chi Squared Distribution
    vary4 = randn(PopulationCases,1)*.1; %Need for random distribution with .1 Standard Dev
    
    % Initializes variables and starts clock for
    tic
    clear out paramtemp LargestProtest wins
    
    % Run Simulations
    for ii =  1:ncases
        paramtemp = ParamCases{ii}; % Selects Case
        
        % **** This is where the case is modified as needed ****
        vary2val = vary2*paramtemp.NVVictoryPercentSTD + paramtemp.NVVictoryPercentMean;
        vary3val = vary3*paramtemp.DefectThresholdSTD + paramtemp.DefectThresholdMean;
        vary4val = round(vary4+paramtemp.PillarProxStrategyMean);
        vary4val(vary4val>1) = 1;
        vary4val(vary4val<0) = 0;
        
        % Runs Each Case Multiple Times
        % parfor jj = 1:length(vary) % Uses Parallel Computing Toolbox for speed
        for jj = 1:length(vary) % Does not require Parallel Computing Toolbox
            temp = paramtemp;
            
            % **** This is where the unique run values are selected ****
            temp.PercentFillActivists = vary(jj);
%             temp.NonviolentSuccessPercent = vary2val(jj);
%             temp.DefectThreshold = vary3val(jj);
%             temp.PillarProxStrategy = vary4val(jj);
            
            % Runs Model
            out = ResistanceABMFinal(sim,temp);
            
            % Gets Needed Results from Each Model
            LargestProtest(ii,jj)=out.ActualMaxInProtestPercentage
            wins(ii,jj)=out.victory;
        end
        casecounter = ii % Output to the screen when case complete
    end
    toc % ends clock
    eval(['save EvoAlg',num2str(LoopCounter),'.mat']) % saves data
    
    %% Post Process
    
    % Determine Good Data
    clear GoodData
    if OnlyGoodData
        GoodData.range = (max(LargestProtest,[],2)>2.5);
        GoodData.StepSize = (max(diff(sort(LargestProtest,2),1,2),[],2)<.6);
        GoodData.NoZeros = sum((wins==0),2)==0;
        % Determines if there is more than one step from losses to wins
        for ii = 1:ncases
            [blah ind] = sort(LargestProtest(ii,:),2);
            tempwins = wins(ii,ind);
            GoodData.switches(ii) = sum(abs(diff(tempwins)))>1;
        end
        GoodData.switches = GoodData.switches';
        GoodData.total = GoodData.switches; % Currently only determines Good Data by switching
    else
        GoodData.total = true(1,ncases);
    end
    ind2plot = find(GoodData.total ==1); %Find Indicies of Good Data
    
    % Make Plots of wins for Good Data
    if PlotRaw
        for ii = 1: length(nonzeros(GoodData.total))
            [blah ind] = sort(LargestProtest(ind2plot(ii),:),2);
            tempwins = wins(ind2plot(ii),:);
            figure
            plot(blah,tempwins(ind),'-x')
        end
    end
    
    % Gets Probability curves using Logistic Regression or Rolling Average
    % Method
    ProtestVector = [];
    WinsVector = [];
    %nSize = [];
    CIlow = [];
    CIhigh = [];
    xpoints = [-1:.5:1.5];
    plusminus = .25;
    nsize = [];
    for ii = 1:ncases
        CompletedRuns = wins(ii,:)>-1;
        if sum(CompletedRuns)>0
            if LogisticOrNot
                xpoints = [-1:.1:1.5];
                temp = GetLogisticRegression(LargestProtest(ii,CompletedRuns)',[wins(ii,CompletedRuns)==2]',xpoints,LogOrNot,0);
            else
                xpoints = [-1:.5:1.5];
                temp = GetRollingAverage(LargestProtest(ii,CompletedRuns),wins(ii,CompletedRuns)==2,xpoints,LogOrNot,plusminus,0);
                %nSize(ii,:) = temp.nsize;
            end
        else
            xpoints = [-1:.1:1.5];
            temp.Y = xpoints.*nan;
            temp.CIlow = temp.Y;
            temp.CIhigh = temp.Y;
        end
        temp2 = GetRollingAverage(LargestProtest(ii,CompletedRuns),wins(ii,CompletedRuns)==2,PopulationBins,LogOrNot,.25,LogOrNot);
        WinsVector(ii,:) = temp.Y;
        CIlow(ii,:) = temp.CIlow;
        CIhigh(ii,:) = temp.CIhigh;
        ProtestVector(ii,:) = xpoints;
        nsize(ii,:) = temp2.nsize;
    end
    
    % Make Population Graph
    figure
    HistoricalData.PopScaled = HistoricalData.BinNumber./sum(HistoricalData.BinNumber);
    PopScaled = nsize./repmat(sum(nsize,2),1,size(nsize,2));
    plot(PopulationBins',PopScaled(GoodData.total,:)','-x',PopulationBins,HistoricalData.PopScaled,'-o')
    ylim([0 1])
    xlim([-1 1.5])
    %xlim([0 3])
    legend('location','SouthEast')
    
    % Make Success Graph
    figure
    plot(ProtestVector(GoodData.total,:)',WinsVector(GoodData.total,:)','-x',HistoricalData.X,HistoricalData.Y,'-o')
    ylim([0 100])
    xlim([-1 1.5])
    %xlim([0 3])
    legend('location','SouthEast')
    
    % Find Best Case for Victory
    MatchUpToVector= xpoints<=MatchUpTo;
    InterpHistoricalWins = interp1(HistoricalData.X,HistoricalData.Y,xpoints,'pchip','extrap');
    error = sum((WinsVector(:,MatchUpToVector)-repmat(InterpHistoricalWins(MatchUpToVector),ncases,1)).^2,2);
    SStot = sum((InterpHistoricalWins(MatchUpToVector)-mean(InterpHistoricalWins(MatchUpToVector))).^2);
    R2 = 1- error/SStot
    
    ErrorPop = sum((PopScaled-repmat(HistoricalData.PopScaled,ncases,1)).^2,2);
    SStotPop = sum((HistoricalData.PopScaled-mean(HistoricalData.PopScaled)).^2);
    R2Pop = 1- ErrorPop/SStotPop
    
    %plot(ShortProtestVector(GoodData.total&(R2Pop>-2),:)',[nsize(GoodData.total&(R2Pop>-2),:)]'/PopulationCases,'-x',HistoricalData.Xfornsize(5:end),HistoricalData.PopScaled,'-o')
    
    TotalError = error/SStot + ErrorPop/SStotPop;
    %TotalError = ErrorPop/SStotPop;
    
    if sum(R2Pop>=R2PopGoal)>0
        ErrorTemp = error;
        ErrorTemp((GoodData.total==0) | R2Pop<R2PopGoal)=nan;
        [MinError MinIndex] = min(ErrorTemp);
    else
        [MinError MinIndexISH] = min(ErrorPop(GoodData.total,:));
        MinIndex = ind2plot(MinIndexISH);
    end
    
    %LowestR2 = min([R2 R2Pop],[],2);
    %[MinError MinIndexISH] = max(LowestR2(GoodData.total,:));
    %MinIndex = ind2plot(MinIndexISH);
    
    % Print and Graph Best Case
    BestParam = ParamCases{MinIndex}
    BestR2 = R2(MinIndex)
    BestR2Pop = R2Pop(MinIndex)
    figure
    handles = plot(ProtestVector(MinIndex,:)',WinsVector(MinIndex,:)','-x',ProtestVector(MinIndex,:)',CIlow(MinIndex,:)','-g',ProtestVector(MinIndex,:)',CIhigh(MinIndex,:)','-g',HistoricalData.X,HistoricalData.Y,'-o',HistoricalData.X,HistoricalData.CIlow,'-c',HistoricalData.X,HistoricalData.CIhigh,'-c');
    ylim([0 100])
    xlim([-1 1.5])
    %xlim([0 3])
    legend([handles(1) handles(2) handles(4) handles(5)],'Best Fit Data','Best Fit CI','Historical Data','Historical CI')
    legend('location','southeast')
    figure
    plot(PopulationBins,PopScaled(MinIndex,:)','-x',PopulationBins,HistoricalData.PopScaled,'-o')
    ylim([0 1])
    xlim([-1 1.5])
    %xlim([0 3])
    legend('Model Data','Historic Data','location','NorthEast')
    
    
    %% Finish Loop
    % Create Next Generation
    if LoopCounter == 1
        StdDevFactor = 2;
    else
        StdDevFactor = 3;
    end
    NewParamCases = newparam4evolution(BestParam,ncases,p,StdDevFactor);
    
    % Loop Variables
    LoopCounter = LoopCounter+1;
    LoopCondition = (LoopCounter <= MaxIter) &  ((BestR2Pop < R2PopGoal) |  (BestR2 < R2Goal));

end